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A lot of effort in lattice simulations over the last years has been devoted to studies of the QCD 
deconfinement transition. Most state-of-the-art simulations use rooted staggered fermions, while 
Wilson fermions are affected by large systematic uncertainties, such as coarse lattices or heavy 
sea quarks. Here we report on an ongoing study of the transition, using two degenerate flavours 
of nonperturbatively 0(a) improved Wilson fermions. We start with N t = 12 and 16 lattices and 
pion masses of 600 to 450 MeV, aiming at chiral and continuum limits with light quarks. 
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1. Introduction 

The high temperature transition from hadronic matter, surrounding us in the universe, to the 
deconfined state of the quark gluon plasma plays an important role in the current state of research 
in nuclear and particle physics, from the theoretical as well as from the experimental point of view. 
Since a lot of experimental effort is being invested in heavy ion collision programmes (CERN SPS, 
LHC and RHIC), investigating the properties of the phase transition and the deconfined phase by 
theoretical means becomes mandatory. This particularly concerns the order of the phase tr ansition 
and the role of the related chiral symmetry restoration. Since the quark gluon plasma is still a 
strongly coupled system for intermediate temperatures, a perturbative treatment of the plasma in 
that region is not valid and lattice QCD is the preferred tool to study the transition. 

Most of the simulations so far were performed using staggered fermions, having the advantage 
of being numerically cheap compared to other fermion discretisations. Recent results can be found 
in refs. [1,2]. However, in order to realise the physical number of quark flavours, fractional powers 
of the staggered quark determinant have to be taken (the "rooting" procedure). The long-standing 
controversy about the validity of rooting (see e.g. [3]) is on-going. Furthermore, the definition of 
the physical pion mass in the rooted staggered approach is ambiguous: At this conference, it was 
debated whether it is correct to refer to the lightest pion mass determined in simulations as the phys- 
ical pion mass, or whether one should use the averaged pion mass as an estimate [4]. This demon- 
strates that a cross-check of the staggered results is needed using other fermionic discretisations. 
Several groups have already started to perform simulations with two flavours of non-perturbatively 
0(d) improved Wilson fermions, of the Sheikholeslami-Wohlert type [5, 6] as well as with the 
maximally twisted mass type [7]. All these simulations still suffer from uncontrolled systematic 
effects, namely unphysically large pion masses and lacking continuum extrapolations, even though 
recent results significantly improved the situation [8, 9]. 

In this paper we report on an ongoing study of the Nf = 2 phase transition with non-pertur- 
batively 6 (a) -improved Wilson fermions, starting at N t = 12 and pion masses m n < 600 MeV. The 
aim is to extract the transition temperature and the order of the transition in the chiral limit, which 
is still not settled until today. There are two possible scenarios: In the first scenario, the chiral 
critical line in the {m U ( t,m s , r}-parameter space never reaches the m Uy d = axis, while the second 
implies the existence of a tricritical point at m u ^ = 0. This tricritical point has to extend into the 
direction of finite chemical potential as a critical line. Thus our Nf = 2 simulation addresses a 
question on the phase diagram of the Nf = 2 + 1 theory and is important for the phase diagram at 
finite density as well. First results suggest that the second scenario is favored by lattice data [5,8], 
but systematics have to be improved in order to make a precise statements. 

2. Setup of the simulations 

2.1 Lattice action and update algorithm 

In our simulations, we work with two degenerate flavours of nonperturbatively &{a) improved 
Wilson fermions, using the Sheikholeslami-Wohlert lattice Dirac operator [10] 

i ci ic * 

D sw = D W + csw —7- cfyy F^ v = D W + csw Msw ■ (2. 1) 
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Here D\y is the usual Wilson Dirac operator, K is the hopping parameter, o^ v the totally antisym- 
metric tensor and F^ v the 'clover leaf representation of the gluonic field strength tensor on the 
lattice. To ensure 0{a) improvement for spectral quantities, the clover coefficient csw has to be 
tuned properly with /3, using the well known interpolation formula for 5.2 < /3 < 12.0 from [1 1]. 

For simulations with large lattices and small pion masses, an efficient algorithm is needed. 
Here we use the deflation accelerated DDHMC algorithm, introduced by Luscher [12]. This al- 
gorithm is nowadays widely used in large-scale lattice simulations and known to work well for 
large lattices up to 128 x 64 3 and pion masses down to 250 MeV (see e.g. [13, 14]). Although the 
algorithm performs well for a variety of lattice sizes it introduces some limitations for simulations 
at non-zero temperature. In particular, we are restricted to lattices sizes N t)S = 8, 12, 16,20,24, . . . 
in each direction with the current version of the algorithm, due to the division of the lattice into an 
even number of DDHMC domains. Following these geometrical constraints and the fact that csw 
is known nonperturbatively only for /3 > 5.20, we see that we have to go to lattices with N t > 12, 
in order to simulate at small pion masses. 

2.2 Observables and strategy 

We investigate the behaviour of three observables around the critical temperature, the average 
plaquette (P), the real part of the Polyakov loop Re [(L)] and the chiral condensate (yw)- In [15] it 
was found that observables constructed from smeared gauge links can lead to a more pronounced 
signal in investigations of the phase structure. Motivated by this reference we also computed the 
Polyakov loop using 5 levels of APE smearing with a weight of 0.5 for the staples. We define the 
generalised susceptibilities %{0) by 



where O is any of the observables above and N s the spatial lattice size. These generalised suscep- 
tibilities should show a notable peak at the transition point. In addition, the behaviour of the peak 
under a change in the spatial volume is governed by the corresponding critical exponents, giving 
information about the order of the transition. 

The temperature T of the system is related to the temporal extent via T = l/(N t a), where a 
is the lattice spacing. There are in general two possible procedures to perform a temperature scan. 
The first is to fix the lattice spacing a and vary the temporal extent N t of the lattice. This procedure 
has the advantage of leaving the scale invariant, but the clear disadvantage is that the resolution is 
limited. This is particularly crucial around the critical temperature and even made worse by the 
use of DDHMC, because of the restrictions on N t discussed above. The other method which we 
use here is to vary the lattice spacing a, while keeping the temporal extent N t fixed. One scheme 
to locate the transition is used by the QCDSF-DIK collaboration: they cross the transition point by 
varying the hopping parameter K, which is related to the bare quark mass, while keeping the lattice 
coupling /3 fixed [5]. The scheme we use here is to fix K while varying /3, enabling us to get a fine 
resolution around T c . Furthermore we are able to employ the multi-histogram method discussed in 
the next subsection. 




(2.2) 
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2.3 Multi-Histogram methods with an explicitly /3 dependent fermion determinant 

Multi-Histogram methods are a tool known for quite some time that allows to use a set of 
measurements at some values {/3 fl } for an interpolation of the expectation values of observables at 
some different value /3 [16]. The method is applicable to theories including fermions, as long as the 
Dirac operator does not depend explicitly on the coupling j3 . For unimproved Wilson fermions this 
is true while for &{a) improved Wilson fermions the determinant becomes jS-dependent. For the 
Sheikholeslami-Wohlert Dirac operator (2.1), the /3-dependence is encoded in the tuning of csw, 
spoiling the applicability of the usual method. In a future publication [17] we will show in detail 
how the multi-histogram method can be extended to the case of &{a) improved Wilson fermions. 
Here we only present an abridged discussion of the main ideas. 

For a /3 -dependent Dirac operator and two degenerate flavours, the weight of a given configu- 
ration [U], generated at /3o £ {/3 fl } for a measurement at /3 is given by 

P(fi)[U] = e-^^- s ^°^R{l5,po)[U}P{l5o)[U] , (2.3) 

where P(j3o) [U] is the statistical weight of the configuration in the j8o-ensemble and we have intro- 
duced an abbreviation for the ratios of fermion determinants at different /3's: 

To have an efficient method, it is necessary to evaluate these ratios without the need for full knowl- 
edge of the determinant, which would demand a new calculation for each value of /3. Using 
eq. (2.1), we can rewrite DswiP ) as 

Dsw(P)[U]=Dsw(fc)[U]+tesw(P,Po)Msw[U], (2.5) 

with Acsw{P,Po) = csw(P) — cswiPo)- Inserting this expression into eq. (2.4) leads to 

R(fi,Po)[U] = (det(l+Ac w (j3, j 8 )D 5V { / (A))[t/]M 5W [[/])) 2 . (2.6) 

In typical /3-scans Acsw <^(10~ 2 ) over the whole range of /3, which suggests that it may be 
sufficient to expand eq. (2.6) in Acsw to approximate the ratio of the two determinants. A short 
calculation yields 

R(fi,Po)[U] = \+2Ac sw (l5,l5 o )Tr(D s ^(p Q )M sw [U]) + 0(Ac 2 sw ) 

= 1+2 Acsw 05 , Po) Rim (ft) [U] + 0{Ac 2 sw ) . 

To apply the multi-histogram method to the full set {A/}> one has to introduce a global density of 
states as in [16]. This is not possible without full knowledge about the determinant at each j5 a . A 
possibility to circumvent this problem, is to introduce a reference point /3r for the fermionic part, 
together with the corresponding global density ^(J3r). In this way, the multi-histogram method 
can be employed using the linearised expression for /?(/3,j3o) discussed above, provided that the 
coefficients multiplying the higher-order terms are not larger than one. It is then sufficient to 
compute the quantity Ru„ (/3 ) for each j8 fl -value of the scan. The combination D^Msw which enters 
Run can be determined with relatively little computational overhead using stochastic estimators. 
Furthermore, Ru n , being the derivative of 7?(j8,j8o) at coupling j8o, provides some control over the 
quality of the linearised ansatz for the ratio of determinants. 
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scan 


Lattice 


Block size 


K 


/3 -range 


T 


tint [f] fmeas 


Statistics 


A 


12 x 24 3 


6 4 


0.13595 


5.27-5.32 


2.0 


0(30) 1 


^(25000) 


B 


16x32 3 


8 4 


0.13650 


5.40-5.55 


2.0 


0(10) 2 


^(4000) 



Table 1: Run parameters for scans A and B. We show the DDHMC block size, the Monte Carlo time T of the 
trajectories, the measurement frequency fmeas an d the integrated autocorrelation time T,,„ of the plaquette P. 



3. Simulation results 

So far our simulations where done on two different lattices, namely 12 x 24 3 (which we refer 
to as scan A) and 16 x 32 3 (scan B). The corresponding parameters are given in table 1. 

Scan A was performed by setting the hopping parameter K to the value where the transition 
in the real part of the Polyakov loop was observed for jS = 5.29 in [5]. On the top of figure 1 we 
show the behaviour of Re [(L)] for scan A. The data are consistent with a phase transition at j8 c = 
5.301(3), where we see a strong increase in the signal and a peak in the susceptibility. By direct 
comparison, we see that the resulting transition temperature is slightly higher than the one obtained 
in [5] . A similar observation was made when comparing the physical transition temperatures from 
[5] to a simulation using twisted mass fermions [9]. Since for the relatively large pion mass of 
roughly 600 MeV the transition is hardly a shaip phase transition but a broad crossover, it is to 
be expected that estimates for the transition temperature derived from different observables or 
directions in parameter space differ. Since we did not perform any scale setting runs at T = 0, we 
do not give an estimate for the transition temperature in physical units. 

Conceptually scan A only served as a feasibility test for our algorithm and our measurement 
routines and an easy tuning case for the simulations at N t > 16. In addition it turned out that this 
scan is at the lower boundary of the working region of the algorithm in the sense that the small 
lattice size limits the possibilities to tune the DDHMC parameters. One of the consequences are 
the relatively large integrated autocorrelation times for the plaquette as can be seen in table 1 , which 
are also reproduced by other observables. These large autocorrelation times are most likely due to 
the relatively small DDHMC blocks of 6 4 , limiting the number of active links during the HMC 
evolution (see [12] for further explanations). Furthermore the pion mass of scan A is still very 
large. We aim at much smaller masses where we also expect a stronger signal for the transition. 

Scan B is our first scan at lighter pion mass and our first attempt for simulations at N t = 16. We 
show the behaviour of Re [(L)] together with the smeared version Re [(L) sm ], and the susceptibility 
of the latter, together with a Gaussian fit to 3 points around the peak position, at the bottom of 
figure 1. The peak position is obtained from the fit as p c = 5.499(2). Fortunately there already 
exists a run for T = with the parameters /3 = 5.50 and K = 0.13650 [13], leading to a transition 
temperature in physical units of roughly T c (m K = 457 MeV, a = 0.060 fm) « 206 MeV. Since the 
determination of the scale is still in progress, we do not give an estimate for the error, but it should 
be kept in mind that the systematical error might still be large. It is important to note that the peak 
is reproduced by all other observables as well. 

4. Conclusions and outlook 

In this proceedings article, we gave a first overview about our ongoing study of the QCD 
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Figure 1: TOP: Plot of the results for scan A. Left: Real part of the Polyakov loop. Right: Correspond- 
ing susceptibility. BOTTOM: Plot of the results for scan B. Left: Real part of the Polyakov loop in the 
unsmeared and APE-smeared version. The smeared results are rescaled such that we could show them in a 
single plot. Right: The susceptibility of the APE-smeared Polyakov loop, together with a gaussian fit at the 
peak position. 



deconfinement transition. We introduced the setup and the main observables to extract the transition 
temperature and an extension of the multi-histogram method to increase the resolution around T c . 
We discussed first results, including the test scan at N t = 12 and a first scan at A 7 , = 16. The 
transition temperature for the latter is at around 206 MeV, but at a still quite large pion mass of 
roughly 450 MeV. The peak in the susceptibilities is reproduced by all observables. The N t = 16 
scan is still extended at the moment to improve the signal and a final analysis is postponed to the 
next publication. We are going to enlarge our set of scans at N t = 16 and also aim at even larger 
values of N t . After mapping the critical trajectories fi c (K,N t ), we are going to do simulations on 
different volumes in order to perform a finite size scaling analysis. In the long run, we aim to do an 
extrapolation to the chiral limit and a separated study of the continuum limit supported by several 
temporal lattice extents and several pion masses for each N t . 
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